#include "field.h"

void field_laplacian_rhu(Field        *f_laplacian,
						 const Field  *f_source,
						 const Field  *f_rho, 
						 const double  inv_dx2,
						 const double  inv_dy2
						 )
{
    
    CHECK_FIELD(f_laplacian);
    CHECK_FIELD(f_source);
    CHECK_FIELDS(f_laplacian,f_source);
    {
        const int       nx        = f_laplacian->nx;
        const int       ny        = f_laplacian->ny;
     //   const int       nx1       = nx-1;
      //  const int       ny1       = ny-1;
	//	const int       nx2       = nx-2; 
        double **       laplacian = ARRAY_OF(f_laplacian);
        const double ** f         = CONST_ARRAY_OF(f_source);
		const double **rho        = CONST_ARRAY_OF( f_rho); 
        int i,j;
		int jp, jm, im, ip, ipp;
        
        for(j=0;j<ny;j++)
		{
			jp=j+1;
			jm=j-1;
			if(j==0)
			{
				jm=ny-1;
			}
			if((ny-1)==j)
			{
				jp=0; 
			}
			for(i=0;i<nx;++i)
			{
				ip=i+1;
				ipp=i+2; 
				im=i-1;
				if(i==0)
				{
					im=nx-1;
				}
				if(i==nx-2)
				{
					ipp=0; 
				}
				if(nx-1==i)
				{
					ip=0;
					ipp=1; 
				}
				const double mid = -2.0 * f[j][i];
				laplacian[j][i]=(
								 inv_dx2 * ( f[j][ip]/(0.5*(rho[j][ipp]+rho[j][ip])) + mid/(0.5*(rho[j][ip]+rho[j][i])) + f[j][im]/(0.5*(rho[j][i]+rho[j][im])) ) +
								 inv_dy2 * ( f[jp][i]/(0.5*(rho[jp][i]+rho[jp][ip])) + mid/(0.5*(rho[j][ip]+rho[j][i])) + f[jm][i]/(0.5*(rho[jm][i]+rho[jm][ip])) ));
			}
		}
    }
}
